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Abstract 

Brownian motion of single particles with various masses M and diameters D is studied by molecular dy- 
namics simulations. Besides the momentum auto-correlation function of the Brownian particle the memory 
function and the fluctuating force which enter the generalized Langevin equation of the Brownian particle 
are determined and their dependence on mass and diameter are investigated for two different fluid densities. 
Deviations of the fluctuating force distribution from a Gaussian form are observed for small particle diame- 
ters. For heavy particles the deviations of the fluctuating force from the total force acting on the Brownian 
particle decrease linearly with the mass ratio m/M where m denotes the mass of a fluid particle. 
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1. Introduction 

Brownian motion has a long history probably 
going back to the end of the eighteenth century 
when Jan Ingcn-Housz 1 mentioned the observation 
of small particles immersed in a fluid that are in 
motion as if they were living beings although not 
the slightest vital spark was in them [1]. Ingen- 
Housz presumably misinterpreted this first docu- 
mented observation of Brownian motion as a conse- 
quence of a fluid motion caused by the evaporation 
of the considered fluid droplet 2 . 

In the first half of the nineteenth century Robert 
Brown published his investigations in which he 
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X J. Ingen-Housz was a dutch scholar and famous physi- 
cian of his time who vaccinated the royal families of king 
George III and of empress Maria Theresia against smallpox. 
He was the first to recognize the importance of chlorophyll, 
the green matter (griine Materie), as he called it, for the 
production of oxygen by plants. 

2 He writes: Wenn man sich auch begniigen wollte, die 
Gcstalt und Grosse von einigen dieser Korpcr wahrend 
des kurzen Zeitraums, als das dem Brennpunkte eines 
Vergrosserungsglases ausgesetzte Tropfchen daucrt, zu 
bcobachten, so mufi man dennoch cingestehen, daB wahrend 
der ganzen Zeit, als das Tropchen dauert, dessen bestandige 
Verdiinstung nothwendigcrweise den ganzen Saft, und fol- 
Preprint submitted to Elsevier 



thoroughly demonstrated that the observed motion 
of small immersed particles is not restricted to liv- 
ing objects but also occurs for inanimate objects as 
long as they are small enough [3]. 

The theoretical understanding of Brownian mo- 
tion by Sutherland [4], Einstein [5] and Smolu- 
chowski [6] opened several completely new vistas 



glich auch die darin cnthaltcnen Korper in eine im- 
merwahrende Bcwcgung versetze, und dafi diese Bewcgung 
betriigen, und in einigen Fallen gewisse Korperchen als 
lebendige Wesen darstellen konne, die nicht den geringsten 
Funkcn des Lebens haben. Um es klar einzusehen, daB man 
sich aus Mangel der Aufmerksamkcit in seinem Urtheile hier- 
iiber betriigen konnte, darf man nur in den Brennpunkt eines 
Mikroskops einen Tropfcn Wcingcist sammt ctwas gestoBcncr 
Kohle setzen; man wird diese Korperchen in einer verwirrten 
bestandigen und heftigen Bewegung erblicken, als wenn es 
Thicrchen warcn, die sich reissend untereinander fortbewe- 
gen. [1, 2]. A very free translation into English reads: The 
study of small particles (literally translated: bodies) im- 
mersed in a small droplet under a magnifying glass is ham- 
pered by two effects: The observation time is seriously lim- 
ited by the evaporation of the droplet (it was Ingcn-Housz 
who later in this paper suggests the use of thin cover slips 
that not yet existed at his time); moreover, the evaporation 
leads to a fluid flow within the droplet that also lets con- 
stantly move the immersed particles such that they seem to 
be alive although not the slightest vital spark is in them. 
This illusion becomes most obvious if one adds under a mi- 
croscope some ground coal to a drop of alcohol; one will then 
see particles in a rapid and permanent motion as if they were 
small animals that violently moved each other. 
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[7]. Brownian motion has provided a cornerstone 
of the mathematical theory of stochastic processes 
in general and of diffusion processes in particular 
[8]. In physics Brownian motion has been a con- 
tinuous inspiration for the understanding of such 
diverse phenomena as transport processes in con- 
densed matter [9], statistical mechanics of non- 
equilibrium processes [10, 11, 12], activated rate 
processes [13, 14, 15], stochastic resonance [16], 
Brownian motors [17, 18], as well as in nanosciences 
[19] and biophysics [20], to name but a few. 

In particular, the description of Brownian motion 
by Langcvin [21] in terms of a Newtonian equation 
of motion including damping and random forces 
has provided an extremely powerful tool to model 
and also to numerically simulate such processes ever 
since we have powerful computers at hand. 

Mori suggested microscopic derivations of so- 
called generalized Langevin equations [22] that also 
allow for memory effects. The Mori formalism 
is based on projection operator techniques that 
can also be used to derive generalized nonlinear 
Langcvin equations [23] as well as corresponding 
generalized master equations [24, 25]. These pro- 
jection operator techniques generally provide in- 
sight into the general structure of the governing 
mesoscopic equations and their symmetries [26] in a 
transparent way but give expressions for transport 
coefficients, memory functions, fluctuating forces 
and the like, that often are extremely involved and 
difficult if not impossible to evaluate apart from cer- 
tain limiting situations [27]. In the particular case 
of Brownian motion of a heavy particle of mass M 
interacting with a large number of fluid particles of 
mass m a Markovian description of the motion of 
the heavy particle results if the mass ratio m/M 
tends to zero [28]. 

On the other hand, Molecular Dynamics (MD) 
provides a direct tool to determine the motion of a 
Brownian particle in a fluid by numerically solving 
the Hamiltonian equations of motion for the Brow- 
nian and fluid particles all of which may interact 
with each other [29, 30]. In this way the friction 
coefficients were determined for an infinitely mas- 
sive particle immersed in truncated Lcnnard- Jones 
and hard sphere fluids in Refs. [31] and [32], respec- 
tively. The assumption of infinite Brownian particle 
mass considerably simplifies the analysis because 
then the velocity of the Brownian particle vanishes 
and the total force acting on the Brownian parti- 
cle coincides with the fluctuating force. Therefore 
the fluctuating force entering the Langevin equation 



becomes directly accessible in the MD simulation 
when the Brownian particle mass is infinite. 

In the present paper we study a single Brown- 
ian particle of various finite masses and diameters 
immersed in a truncated Lennard- Jones fluid. The 
range of parameters that we consider here is rel- 
evant for diffusion of medium size molecules like 
fullcrenes in liquids. For applications to nanoflu- 
idics the influence of the confinement becomes im- 
portant but is not considered here. 

In the studied parameter regime the fluctuating 
force in general differs from the total force by the 
so-called organized time rate of change of the mo- 
mentum which contains a memory function weight- 
ing the influence of the value of the momentum as it 
was at earlier times. This memory function though 
cannot directly be determined from the output of 
the MD simulations which consists of the position 
and momentum of the Brownian particle, as well as 
of the total force acting on the Brownian particle, 
at each integration time step. In fact, the memory 
kernel is related to the momentum auto-correlation 
function by a Volterra integral equation of the first 
kind [10]. The numerical solution of this type of 
equations is known to be error-prone. A conversion 
to a Volterra equation of the second kind can be 
achieved, for which more stable algorithms exists, 
but which involves derivatives of the momentum 
auto-correlation function with respect to time up 
to second order. For their efficient estimates Berne 
and Harp [33] used local polynomial approxima- 
tions of the momentum auto-correlation function. 
Here we do without any fitting. We identify the 
first and second time derivatives of the momentum 
auto-correlation function with correlation functions 
of momentum and total force and with the total 
force auto-correlation function, respectively. Both 
correlation functions can be directly estimated from 
the MD simulation data leading to reliable results 
for the memory function. 

Kncller and Hinsen suggested an alternative 
method of determining the memory function [34]. 
Their method is based on a fit of the momen- 
tum auto-correlation function by an auto-rcgrcssive 
model AR(P) where the order P determines the 
maximal extent of the memory time. We do not 
further pursue this method since it requires the fit- 
ting of a large number of auxiliary parameters. 

Once the memory kernel is known, the fluctuat- 
ing force can also be determined and further ana- 
lyzed. Of particular interest is the questions under 
which conditions the random force becomes Gaus- 



sian and/or Markovian and when it can be approx- 
imated by the total force. 

The paper is organized as follows. In Section 2 
the employed microscopic model is specified and the 
MD simulation outlined. In Section 3 the general- 
ized Langevin equation is reviewed and the esti- 
mation of the respective memory kernel and fluc- 
tuating force based on MD simulation data is de- 
scribed. The momentum auto-correlation function 
obtained from the MD simulations is discussed in 
Section 4. The consistency of the resulting mem- 
ory kernel with the structure of the generalized 
Langevin equation imposed by the Mori theory is 
confirmed in Section 5. The mass and diameter de- 
pendence of the memory kernel and the fluctuating 
force is investigated in Sections 6 and 7. The paper 
ends with concluding remarks. 

2. Microscopic model of 2d Brownian mo- 
tion and MD simulations 

As a simple microscopic model of Brownian mo- 
tion we consider N soft "fluid" particles of mass m 
and diameter d and a single "Brownian" particle 
of mass M and diameter D moving in a two di- 
mensional quadratic domain of side-length L with 
periodic boundary conditions. With m = M and 
d = D this model includes the case of self-diffusion. 
The typical regime of Brownian motion though is 
described by M > m and D 3> d. The fluid par- 
ticles interact pairwise with each other as well as 
with the Brownian particle. The Hamiltonian de- 
scribing the classical motion of this N + 1 particle 
system is of the form 



1 N 1 N 
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where and Q denote the positions of the fluid and 
Brownian particles, respectively, and Pi and P the 
according momenta. The interaction V(r) is purely 
repulsive and equally acts between pairs of fluid 
particles and between the Brownian particle and 
fluid particles. It is given by a truncated, purely 
repulsive Lcnnard-Jones potential of the form 



V<r(r) 



(cr/r) -O/r) 



e for r < 2 1 / 6 a 
for r > 2 1 / 6 cr 
(2) 



In the MD simulations dimensionlcss units were 
used for which the fluid particle assume the diame- 
ter d = (Tfi = 1. For Brownian particles of diameter 
D we chose <jb = {D + d)/2. Hence, <jb is the con- 
tact distance between the Brownian particle with 
diameter D and a fluid particle with diameter d. 
Masses are given in multiples of the fluid particle 
mass and the energy unit e is chosen as difference 
of the pair potential energy at distances r = aji 
and r = 2 1 / 6 a // , e = V afl {a fl ) - V ajl {2^ 6 a fl ). A 
consistent unit of time is then given by r = d^Jm/e. 

Then the reduced parameters length r* , temper- 
ature T* , and time t* are defined as r* = r/d, 
T* = Tk B /e, and t* = t/r. If if is defined as 
the ratio of the area occupied by the fluid particles 
to the total available area 77* = Nird 2 /(4L 2 ), the 
reduced density is defined as n* = 4?f /it. We have 
considered a fluid system composed of N = 10, 000 
fluid particles and one Brownian particle. The sim- 
ulations were performed at two different densities, 
one at n* = 0.4, and the other one at n* = 0.8. 
The temperature was always T* = 1. 

Initially the fluid particles occupy lattice points 
of the (111) face of an fee lattice. To each particle 
a random two dimensional vector p° = ^/2mfcBTe; 
with identically, uniformly distributed unit vectors 
e.j is assigned. By subtracting the average p° = 
TV -1 YliLi Pi we generated initial values of the fluid 
particle momenta p^ = p° — p° such that the total 
fluid is at rest. 

The initial position and momentum of the Brow- 
nian particle both are zero. The simulations were 
realized at constant energy using the standard ve- 
locity Verlet algorithm [29] with a time step h = 
10 _3 t to insure the stability of the total energy to 
within 10~ 3 %. 

We found that it never took more than 5 x 10 6 
time steps until thermal equilibrium for fluid and 
Brownian particle was established. From the subse- 
quent 3 x 10 7 time steps every third one was used for 
the reconstruction of a generalized Langevin equa- 
tion as described in the rest of this paper. 



3. Generalized Langevin equation and its es- 
timation from molecular dynamics simu- 
lations 

3.1. Mori's generalized Langevin equation 

The Mori theory provides a framework to deter- 
mine the equilibrium correlation functions of a set 
of so-called "relevant" or "macroscopic" variables 



as well as the relaxation of their mean values close 
to equilibrium in terms of generalized Langevin 
equations. These equations express the time rate of 
change of the relevant variables as the sum of sys- 
tematic and random contributions, the latter also 
known as the random force. The relevant part is 
a linear expression in the relevant variables and in 
general consists in an instantaneous reversible and a 
retarded contribution. The retarded contribution is 
determined by a memory kernel that is connected to 
the auto-correlation function of the random part via 
a fluctuation dissipation theorem of second kind. 
The correlation of the random force with the rel- 
evant variables vanishes. Higher order correlation 
functions of the fluctuating force and consequently 
higher than second moments of the relevant vari- 
ables remain unspecified within Mori theory which 
hence does not provide a complete characterization 
of the stochastic process of a given set of relevant 
variables. 

In the present study of free Brownian motion we 
choose the momentum P of the Brownian parti- 
cle as the relevant variable. Due to the absence 
of an external potential there is no instantaneous 
contribution to the systematic part of the momen- 
tum time rate of change. Hence, the generalized 
Langevin equation takes the form 

P(t) = - [ dsk(f-s)P(s)+F + (t), (3) 
Jo 

where the mean value of the fluctuating force van- 
ishes, i.e. 

(F+(t))=0. (4) 

Moreover, the correlation functions of the fluctuat- 
ing force components F a (t) are related to the com- 
ponents of the memory kernel k a ,p(t) by the fluc- 
tuation dissipation theorem 

(F+(t)F+(s)) = J2(P*P^At ~ *) ■ (5) 

7 

Due to the isotropy of the present microscopic 
model specified in the previous section all non- 
diagonal components vanish and the diagonal com- 
ponents agree with each other 

(F+(t)F+(s)) =6 a , f3 {F + (t)F + (s)) , 

(P a Pp) = 6 a ,p(P 2 ) , (6) 
k a ,p(t - s) = 5 a ,pk(t - s) . 

Here F + (t) and P denote either of the x or y com- 
ponents of the fluctuating force and momentum, 



respectively. Consequently, both components sat- 
isfy the same scalar generalized Langevin equations 
reading 

P(t) = -[ dsk(t-s)P(s)+F+(t) , (7) 
Jo 

with F + (t) obeying the scalar fluctuation dissipa- 
tion theorem 

(F+(t)F+(s)) = {P 2 )k(t-s). (8) 

The microscopic derivation of the generalized 
Langevin equation (7) and the corresponding mi- 
croscopic expressions of the fluctuating force and 
the memory kernel are based on an identity for the 
microscopic time evolution operator e~ ct reading 

e -ct = e -ct p + (1 _ P)e -d-v)ct {1 _ v) 

- f ds e- c{t - s) VC(l - V)e- {1 - v)Cs (l - V) , 
Jo 

(9) 

where C denotes the Liouville operator governing 
the time evolution of the microscopic phase space 
probability density p 

dt Lp (10) 

The operator V provides an orthogonal projection 
of phase space functions / onto the linear subspace 
of relevant variables which are the x and y compo- 
nents of the Brownian particle momentum in the 
case of free Brownian motion. It is defined as 

vf = Y,p«(p a j)/(p a ,p a ), (ii) 

a 

where 

(9, f) = J dTgfe- H '^ T / J dTe~ H ^ T (12) 

denotes the Mori scalar product of phase space 
functions g and / with respect to the Maxwell- 
Boltzmann distribution e ~ H / k sT/ J dTe -H/k B T de _ 
scribing thermal equilibrium of fluid and Brownian 
particle at the temperature T. Here the integral 
extends over the phase space with the volume el- 
ement dT = dPdQY\^ dpidqi. Hence, the Mori 
scalar product of two phase space functions agrees 
with the thermal expectation value of the product 
of these functions, i.e. 

(f,g) = (fg)- (is) 



Applying the identity (9) to the time rate of 
change of either component of the momentum given 
by P(t) = —e~ ct CP one obtains the generalized 
Langcvin equation (7) with the following micro- 
scopic expressions for the fluctuating force and for 
the memory kernel: 

F + (t) = (l-T)exp{(l-V)Ct}P , (14) 
k(t) = (Pe- {1 - p)Cs P) , (15) 

where P = —CP. 

By a scalar multiplication of both sides of eq. (7) 
with P(0) one obtains the equation of motion for 
the momentum auto-correlation function 

C(t) = (P(t)P) (16) 

reading 

C{t) = - f dsk(t- s)C(s) . (17) 
Jo 

Thereby one uses the fact that the fluctuating force 
F + (t) and the momentum P(0) are orthogonal, i.e. 

(F+(i),P(0))=0, (18) 

as follows from eq. (14). 

3.2. Memory kernel and fluctuating force from MD 
simulations 

From an MD simulation of a system of fluid par- 
ticles interacting with a single Brownian particle as 
described in Sect. 2, the knowledge of the instanta- 
neous positions and momenta of all particles allows 
one to directly obtain the momentum of the Brow- 
nian particle as well as the total force acting on the 
Brownian particle. However, the separation of the 
total force into an organized and a random contri- 
bution as it is presented in the generalized Langevin 
equation (7) cannot be inferred from the instan- 
taneous microscopic state of the total system. In 
order to achieve this separation one first estimates 
the stationary momentum auto-correlation function 
C(t), see eq. (16) and then uses the equation of 
motion of the momentum auto-correlation function 
(17) in order to determine the memory kernel. 

At first sight it is tempting to employ a Laplace 
transformation which changes the convolution of 
the memory kernel and the momentum auto- 
correlation function into a product and immediately 
leads to an explicit expression for the Laplace trans- 
formed kernel. However, the extreme sensitivity of 



the inverse Laplace transform to numerical errors 
renders a reliable determination of the memory ker- 
nel in the time domain practically impossible. 

Once the momentum auto-correlation function 
is known, eq. (17) represents a Volterra integral 
equation of first kind for the memory kernel. A 
discretization of this equation in principle can be 
solved quite effectively because it involves the inver- 
sion of a Toplitz matrix. However, also this method 
is plagued by numerical inaccuracies. 

More stable algorithms exist for Volterra equa- 
tions of second kind. Any Volterra equation of first 
kind can be transformed to a Volterra equation of 
second kind by differentiation with respect to the 
independent variable. In the present case of eq. 
(17) this yields 

C(t) = -C(0) k(t) - [ dsC(t- s)k(s) . (19) 
Jq 

Here, however, the numerical differentiation of 
the momentum auto-corrclation function, which is 
not analytically known, may introduce large nu- 
merical errors. Berne and Harp [33] attempted 
to keep these errors under control by applying 
fourth and sixth order polynomial approximations 
of the velocity auto-correlation for small and large 
times, respectively [33]. We here choose a differ- 
ent strategy avoiding any numerical differentiation. 
The first and second derivatives of the momentum 
auto-correlation function can be expressed as the 
total-force-momentum correlation function and the 
total-force auto-correlation function, respectively. 
Strictly speaking, we have 

C(t) = (F(t)P) , (20) 
C(t) = -(F(t)F). (21) 

These correlation functions with the total force can 
be directly estimated from the MD simulation and 
therefore do not introduce additional errors. Hence, 
we determined the memory kernel from a numerical 
solution the integral equation 

(F(t)F) = C(0)k(t)+ f ds(F(t-s)P)k(s) , (22) 
Jo 

which follows from eqs.(19),(20) and (21). The nu- 
merical scheme to solve eq.(22) is rather straight- 



forward. The discretization of eq.(22) reads 
(F(iAt)F) 

i 

= C(0) k(iAt) + At^2^j (F{iAt - jAt)P) k{jAt) , 

j'=o 

(23) 

where Wj = 1/2 for j = 0, i and tUj = 1 otherwise 
are a weight factors for the integration. Then, k(t) 
at every t = iAt can be obtained iteratively as 

k(iAt) = {C*(0) + AtLo t {FP)y x x 

{F(i*t)F) -AtJ2 u, (F((i - j)At)P) k(jAt) 1 , 

3=0 ) 

(24) 

with initial condition fc(0) = (F 2 )/C(0). Such ob- 
tained memory kernels for different mass, diameter 
and density values are displayed in the Figs. 6 and 
7 and further discussed in Section 6. 

Once the memory kernel is known, the fluctuat- 
ing force can be calculated as the difference of the 
total force that is taken from the MD simulation 
and the systematic part of the force, i.e we have 

F + (t) = F{t)+ [ dsk(t- s)P(s) , (25) 
Jo 

where the values of the momenta at the times s 
prior to t are also taken from the MD simulation. 
Probability densities of the fluctuating force esti- 
mated from histograms for different masses and di- 
ameters of the Brownian particle and different den- 
sities are displayed in Fig. 8 and further discussed 
in Section 7. 

3.3. Initialization and stationarity 

The Mori equation (7) allows for relaxations of 
the momentum as long as the initial momentum 
is sufficiently small such that it can be described 
within linear response theory. This means that 
the initial probability distribution of the fluid plus 
Brownian particle system must be of the form 

p (Q,P,q.p) = Z -l (Po)e -[H-P.P /*l/* a T 
«Z- 1 (0)e-^ T (l + -lPo-p) , ^ 

where Po denotes the averaged initial momentum 
and Z(Po) the partition function that is denned as 

Z(P ) = f dr e -^- p - Vo ^ hB . (27) 



This corresponds to a situation in which the bath is 
in thermal equilibrium in the presence of a Brow- 
nian particle moving on average with momentum 
Po 

In the MD simulations it though is more practical 
to start with initial conditions in which neither the 
bath nor the Brownian particle are at equilibrium 
and then let run the simulation long enough until 
equilibrium has established. We describe the result- 
ing stationary process by a generalized Langevin 
equation in which the time of the initial prepara- 
tion is shifted to the infinitely remote past, i.e. we 
do not make reference to any particular initial con- 
dition and consider 

P(i) = - f dsk(t - s)P(s) + F+(t) . (28) 

J — oo 

In order to be consistent with the Mori equation 
(7) we require the relation (17) between the mem- 
ory kernel k(t) and the momentum auto-correlation 
to hold. Then the fluctuating force and the momen- 
tum are correlated in the following way 

/>oo 

(F + (t)P)= dsk(t + s)C(s) . (29) 
Jo 

In contrast, in the framework of the Mori equa- 
tion (7) the initial momentum and the fluctuating 
force at positive times are uncorrelated as it fol- 
lows from cq. (18) because of the initial prepa- 
ration (26). For a formal proof of eq. (29) see 
the appendix A. Moreover, one can show that the 
stationary Langevin equation in conjunction with 
the relation (29) implies the fluctuation dissipation 
theorem (5). A proof is also given in appendix A. 
Vice versa, if one imposes the fluctuation dissipa- 
tion theorem and the relation (29) then the equa- 
tion of motion for the momentum auto-correlation 
function, eq. (17), is recovered, as proved in Ap- 
pendix A. Therefore, the Mori equation (7) and the 
generalized Langevin equation (28) provide equiv- 
alent descriptions of the equilibrium properties of 
the Brownian particle. 

4. Momentum auto-correlation functions 

After a transient period, the Brownian parti- 
cle has equilibrated. Accordingly, the momen- 
tum distribution of the Brownian particle becomes 
Maxwellian. In particular, the momentum distri- 
butions are independent of the size of the Brow- 
nian particles and of the fluid density; moreover, 




Figure 1: (color online) The probability density of a momen- 
tum component p{P) = exp{- P 2 / (2Mk B T)}/^2nMk B T 
(solid line) for M = 100m is compared to the histograms 
containing the x- and y-components of the momentum ob- 
tained from MD simulations for Brownian particles with dif- 
ferent diameters and two different fluid densities: D = lOd 
(small (n* = 0.8) and large (n* = 0.4) black plus), D = 7d 
(empty (n* = 0.8) and filled (n* = 0.4) blue diamond), 
D = 5d (small (n* = 0.8) and large (n* = 0.4) red cross), 
D = 2d (open (n* = 0.8) and filled (n* = 0.4) green circles), 
and D = d (open (n* = 0.8) and filled (n* = 0.4) magenta 
square). The agreement between theory and simulation is 
very good. 

the widths of the distributions conform with the 
equipartition law, (P%) = MksT for each compo- 
nent. These findings are illustrated in Fig. 1 for 
Brownian particles of mass M = 100?n and differ- 
ent diameters and densities. The simulation results 
for other masses of the Brownian particle also con- 
form with theory but are not shown. 

The stationary momentum auto-correlation func- 
tion C(t) = (P(t)P(Q)) is estimated as time average 
of the stationary part of the simulated time series 
of momenta, i.e. 

= JM ? P(to+jAt + t)P(t Q +jAt) (30) 
j=i 

where t$ is the time after which equilibrium has 
established. The upper limit of the sum is given 
by Af(t) = N - tj At where TV is the total num- 
ber of momenta. The correlation function was es- 
timated for every third time step At = 3h. Re- 
sults arc displayed for Brownian particles with mass 
equaling the fluid mass and different diameters in 
Figs. 2(c) and 2(d) and for heavy Brownian par- 
ticles in Figs. 2(a) and 2(b). The dependence of 
the normalized momentum auto-correlation func- 
tion C(t)/C(0) on mass is illustrated for small and 
large Brownian particles in Fig. 3, respectively. The 
momentum auto-correlation function of light and 



large particles is characterized by decaying oscilla- 
tions whereas the decay becomes monotonic in the 
case of heavy or small particles. For the lower den- 
sity the period of oscillations is larger and the speed 
of decay is decreased due to a larger mean free path 
and mean time between collisions. The oscillations 
observed for large light Brownian particles indicate 
the presence of a cavity surrounded by fluid parti- 
cles within which the Brownian particle moves back 
and forth. 



An estimate of the large time behavior of the 
momentum auto-correlation function is a notori- 
ously difficult problem [35]. These long-lived cor- 
relations though deteriorate the estimate of the 
momentum statistics as we will discuss now. By 
comparing the distributions of the x and y com- 
ponents of the momentum one finds slight devia- 
tions, see Fig. 4. These deviations can mainly be 
attributed to differences of the estimated mean val- 
ues P a = E?=_iP a (to + a = x,y. The 
time averages P a themselves are random quantities 
having the ensemble average (P a ) = and the vari- 
ancc (P») = (C(0) + 2Y? =1 {N- i)C(iAt))/U 2 « 
2 J Q dt(T — t)C(t) /T 2 . The approximation by an 
integral is valid in the case of oversampled data, 
i.e. if C(t) changes only little with the sampling 
time At. In any case, the variance of the mo- 
mentum time average increases when the correla- 
tions extend over large times, explaining the slight, 
seeming anisotropy of the momentum statistics for 
heavy Brownian particles. For light Brownian par- 
ticles the momentum auto-correlation decays much 
faster, (cf. Figs. 2 and 3), and therefore the esti- 
mated momentum distribution conforms much bet- 
ter with the expected isotropy, see Fig 4. 



A better estimate of the statistics of the momenta 
can be obtained if it is based on several shorter in- 
dependent trajectories instead of a single long one. 
We confirmed this by comparing estimates from 
ten trajectories consisting of 10 6 sampled momenta 
with a single one with 10 sampled momenta. The 
former indeed led to better convergence. However, 
since we here are mainly interested in the mem- 
ory kernel and the statistical properties of the fluc- 
tuating forces, these quantities are essentially de- 
termined by the time correlations of the momenta 
based on single long trajectories. 




Figure 2: (color online) The momentum auto-correlation function C(t) is displayed for Brownian particles of mass M = 100m 
in panels (a) and (b) and M = m in panels (c) and (d). In panels (a) and (c) the fluid density is n* = 0.8, and n* = 0.4 for 
panels (b) and (d). In each panel different diameters D = lOd (red, dashed line), D = 5d (green, solid line), D = 2d (blue, 
dash-dotted line) and D = d (black, dotted line) are shown. All cases shown correspond to the temperature T* = 1. The 
larger density causes more frequent collisions of the Brownian particle with fluid particles and consequently a faster decay of 
the momentum auto-correlation function. The pronounced oscillations of the auto-correlation function of large, light particles 
indicates the presence of an almost empty cavity in which the particle moves back and forth, see panels (c) and (d). The 
period of oscillations is independent of the particle diameter. For smaller particles the damping of the oscillations increases. 
The oscillation period is larger for the low density fluid with n* = 0.4. 
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Figure 3: (color online) The mass dependence of the normalized momentum auto-correlation function C(t)/C(0) is displayed 
for Brownian particles with diameter D = lOd in panels (a) and (b) and for diameter D = d in panels (c) and (d). Panels (a) 
and (c) are for a high fluid density n* = 0.8 and panels (b) and (d) for the lower density n* = 0.4. In all panels the masses 
are M = 100m (red, dashed line), M = 25m (green, solid line), M = 4m (blue, dash-dotted line) and M = m (black, dotted 
line). The pronounced oscillations for large light particles disappear with increasing mass, see (a) and (b), and are absent for 
small particles of radius D = d at any mass, see (c) and (d). Lower density leads to a slower decay and in case of oscillations 
to a longer period. 



5. Two consistency checks 




Q*(P) 



Before discussing the memory kernel and the sta- 
tistical and temporal properties of the fluctuating 
force in more detail we consider the consistency of 
the estimates of these quantities obtained from MD 
simulations by checking the relations (29) and (8). 

The stationary fluctuating force-momentum cor- 
relation function (F + (t)P) entering eq. (29) was 
estimated from MD simulation results for Brown- 
ian particles and compared with the integral of 
the product of the memory kernel shifted in time 
and the momentum auto-correlation function. Also 
these functions were determined from the MD sim- 
ulation data. Fig. 5 displays a perfect agreement 
of the fluctuating force-momentum auto-correlation 
function and this integral for several masses and di- 
ameters of the Brownian particle as well as for two 
densities of the fluid. 

We also find perfect agreement of the fluctuat- 
ing force auto-correlations function normalized by 
C(0) = (P 2 ) with the memory kernels for the same 
set of Brownian particle and fluid parameters as 
above, see Fig. 6. Hence the simulated MD data 
of the Brownian particle momenta are perfectly de- 
scribed by a generalized Langevin equation with the 
specified memory kernel and fluctuating forces re- 
lated by the fluctuation dissipation theorem. 



Figure 4: (color online) The QQ-plot of the y and x mo- 
mentum components displays a slight deviation between the 
cumulative distribution functions Q a (P) = Prob(P Q < P) 
with a = x, y for a heavy Brownian particle with mass 
M = 100m (red, dashed line) in contrast to the case of a 
light particle with M = m (blue crosses). In both cases 
particles have the same diameter D = lOcZ and move in a 
fluid of density n* = 0.8. The ideal relation Q y = Q x indi- 
cating identical distributions is displayed by the thin black 
line. The seeming anisotropy of the momenta is of statisti- 
cal nature due to the long-lived momentum auto-correlation 
function of heavy Brownian particles. 



6. Brownian mass and diameter dependence 
of the memory kernel 

As it is obvious from Fig. 6 the initial value of the 
memory kernel k(0) strongly depends on the mass 
M of the Brownian particle. From the equation of 
motion of the momentum auto-correlation function, 
eq. (17), it follows upon differentiation with respect 
to time that 

(7(0) 



fc(0) 



C(0)' 



(31) 



Assuming that C(0) only weakly depends on mass 
we obtain with the equipartition law C(0) = MksT 
that fc(0) is inversely proportional to the mass of the 
Brownian particle. For a comparison of the mem- 
ory kernel at different masses M we therefore scale 
its value by the mass ratio m/M, i.e., we consider 
the dimensionless reduced memory kernel fc(t)r 2 j 
as a function of where Tyi = Ty/M/m is the 
relevant, "macroscopic" time scale of a Brownian 
particle of mass M. Fig. 7 clearly indicates that 
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Figure 5: (color online) The stationary fluctuating force momentum correlation function {F + (t)P) displayed as lines is compared 
to the integral dsk(t + s)C(s) (symbols) for heavy Brownian particles with Mass M = 100m in panels (a) and (b) and light 
Brownian particles, M = m in panels (c) and (d). Panels (a) and (c) refer to the density n* = 0.8 and panels (b) and (d) to 
the lower density n* = 0.4. All panels contain results for different diameters D = lOcZ (red, dashed line; red circle), D = 5d 
(green, solid line; green cross), D = 2d (blue dash-dotted line; blue diamond) and D = d (black, dotted line; black plus). In all 
cases the agreement is perfect in accordance with eq. (29). 
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Figure 6: (color online) The memory function k(t) displayed as lines is compared to the normalized fluctuating force auto- 
correlation function (F + (t)F+)/C(0) shown by symbols for heavy Brownian particles with Mass M = 100m in panels (a) and 
(b) and light Brownian particles, M = m in panels (c) and (d). Panels (a) and (c) refer to the density n* = 0.8 and panels (b) 
and (d) to the lower density n* = 0.4. All panels contain results for different diameters D = lOd (red, dashed line; red circle), 
D = 5d (green, solid line; green cross), D = 2d (blue dash-dotted line; blue diamond) and D = d (black, dotted line; black 
plus). In all cases the agreement is perfect in accordance with the fluctuation dissipation theorem (8). 
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Figure 7: (color online) The rescaled memory function k(t)T M is displayed in dependence of the macroscopic time t/TM for 
large Brownian particles with diameter D = lOd in panels (a) and (b) and small Brownian particles with D = d in panels (c) 
and (d). Panels (a) and (c) correspond to the fluid density n* = 0.8 and panels (b) and (d) to n* = 0.4. All panels contain 
results for different masses M = 100m (red, dashed line), M = 25m (green, solid line), M = 4m (blue dash-dotted line), and 
M = m (black, dotted line). In the insets the short time behavior of the rescaled memory kernel is shown. Note that the 
rescaled memory function assumes a value at t = that is almost independent of the mass. This confirms the assumption that 
the second derivative of the momentum auto-correlation function only weakly depends on the mass. The decay of the memory 
kernel on the macroscopic time scale though becomes faster with growing mass and also with decreasing diameter, see also 
Fig. 6. 
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the rcscaled memory kernel at zero time is almost 
independent of mass. 

The decay of the memory kernel is qualitatively 
the same for all values of Brownian particle mass 
and diameter as well as for both fluid densities. It 
is much faster than the decay of the correspond- 
ing momentum auto-correlation functions and for 
most parameter values consists in a rapid initial 
decay that overshoots to negative values and then 
slowly approaches zero from the negative side. In 
particular, the memory kernel docs not display pro- 
nounced oscillations in contrast to the momentum 
auto-correlation function. As for the momentum 
auto-correlation function our data are not sufficient 
to specify a particular decay law at large times. 

The overall decay is faster both for heavier and 
smaller particles. Hence, with increasing mass 
of the Brownian particle an approximation of the 
memory kernel by a delta function describing in- 
stantaneous friction becomes more reliable. The 
larger density leads to a larger reduced memory 
function but hardly influences its shape. 

7. Brownian mass and diameter dependence 
of the fluctuating force 

The statistical properties of the fluctuating force 
strongly depend on the diameter of the Brownian 
particles and the fluid density. For small particles 
and low densities the distributions are leptokurtic 
with approximately exponential, or even more pro- 
nounced tails. For larger Brownian particles and 
higher fluid densities the distributions approach a 
more Gaussian shape, see Fig. 8. A distinct depen- 
dence of the fluctuating force distribution on the 
mass of the Brownian particle only exists for small 
Brownian particles and low fluid densities while it 
is insignificant otherwise. 

The observed statistics of the fluctuating forces 
can be qualitatively understood in terms of the 
number of fluid particles that simultaneously in- 
teract with the Brownian particle. This number, 
Ni, fluctuates as fluid particles constantly enter and 
leave the interaction region of the Brownian parti- 
cle. If, for example, at a given instant of time, no 
fluid particle happens to interact with the Brownian 
particle, i.e. if Nj = 0, then there will be no force 
exerted on the particle. This situation may happen 
more frequently if the fluid density is low and the 
Brownian particle diameter and therefore also its 
interaction region is small. This explains the pro- 
nounced peak of the fluctuating force distribution 
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Figure 8: (color online) Histograms of the fluctuating forces 
are displayed for Brownian particles with mass M = 100m 
and different diameters D = d (blue, dotted line) D = 5d 
(red, solid line) and D = lOd (black broken line). In the 
latter case a Gaussian distribution with the same mean value 
and variance is shown for comparison (green, thin solid line). 
It perfectly coincides with the histogram for the high fluid 
density n* = 0.8 in panel (a) whereas deviations are visible 
in the tails of the distribution for the lower density n* = 0.4 
in panel (b). 



at zero force for low fluid density and small Brow- 
nian particles. On the other hand, if the number of 
fluid particles interacting with the Brownian parti- 
cle is large, as it is the case for large densities and 
large Brownian particles, not only the probability 
of zero force decreases but also very large forces re- 
sulting from single particle impacts are likely to be 
almost compensated by the influence of the other 
interacting particles and hence the frequency of ex- 
tremely large forces becomes suppressed. In this 
way probability is transferred both from the tails 
and the center towards the flanks of the force dis- 
tribution when the fluid density or the particle di- 
ameter increase. Since the forces that are simulta- 
neously acted on the Brownian particle by different 
fluid particles will only show weak correlations, the 
sum of these forces will converge to a Gaussian ran- 
dom number if, on average, the number Nj is suffi- 
ciently large. Although these arguments primarily 
apply to the total force F(t) acting on the Brown- 
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Figure 9: (color online) The variance of the fluctuating forces 
is displayed for fluids with densities n* = 0.8 and n* = 0.4 in 
the panels (a) and (b), respectively, as functions of the diam- 
eter of the Brownian particle D for different masses M = m 
(black plus), M = 4m (red cross), M = 25m (green circle), 
M = 49m (blue diamond) and M = 100m (yellow triangle). 
The variances for the different masses differ only insignifi- 
cantly from the behavior obtained for the average over the 
masses (black solid line). The variance of the fluctuating 
forces for the denser fluid is larger, but otherwise the increase 
at small diameters and the saturation to an approximately 
constant large particle value agrees qualitatively. 



ian particle it still reflects the observed behavior of 
the fluctuating force statistics. 

In order to characterize the parameter depen- 
dence of the fluctuating force statistics in more 
quantitative terms we estimated the first four mo- 
ments of the fluctuating force for various parameter 
vales. All odd moments do not significantly differ 
from zero. The second moment therefore coincides 
with the variance of the fluctuating force. It in- 
creases with the diameter D and approaches a con- 
stant value for large diameters. At lower fluid den- 
sities the variance is smaller but it hardly depends 
on the mass M of the Brownian particle. These 
features are illustrated in Fig. 9. 

From the second and fourth moment of the 
fluctuating force the kurtosis is determined as 
(F+ 4 ) / (F+ 2 ) 2 - 3. As it measures the excess of 
probability located in the center and the tails of 
a distribution as compared to a Gaussian distri- 




Figure 10: (color online) The kurtosis (F +i ) / (F+ 2 ) 2 - 3 
of the fluctuating forces is displayed for fluids with densities 
n* = 0.8 and n* = 0.4 in the panels (a) and (b), respectively, 
as functions of the inverse diameter of the Brownian parti- 
cle 1/D for different masses M = m (black plus), M = 4m 
(red cross), M = 25m (green circle), M = 49m (blue dia- 
mond) and M = 100m (yellow triangle). As a guide to the 
eye, points are connected by solid lines of same color. For 
increasing particle size the kurtosis approaches zero. 



bution, it is larger for the less dense fluid and ap- 
proaches zero with increasing diameter of the Brow- 
nian particle, see Fig. 10, in agreement with the 
qualitative arguments discussed above. 

Whereas the kurtosis only specifies a particular 
aspect of a deviation from a Gaussian distribution, 
the Kullback-Lcibler divergence [36] allows one to 
quantify a distance measure of a given distribution 
from the closest Gaussian distribution. According 
to its definition the Kullback-Leibler divergence of 
p* to p is given by 



KL(p\\p*) = J dxp(x)\n-^ 



(32) 



where p and p* denote probability density func- 
tions. The Kullback-Leibler divergence in general 
is positive. It vanishes if and only if the probabil- 
ity density functions p and p* coincide with each 
other [36]. It is straightforward to show that the 
Gaussian distribution p* with the mean value and 
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Figure 11: (color online) Minimal Kullback-Lcibler diver- 
gence of the fluctuating force distribution to a Gaussian dis- 
tribution KL = min Gaussianp . KL(p F+ 1 \p*) is displayed for 
fluids with densities n* = 0.8 and n* = 0.4 in the panels (a) 
and (b), respectively, as functions of the inverse diameter 
of the Brownian particle 1/D for different masses M = m 
(black plus), M = 4m (red cross), M = 25m (green circle), 
M = 49m (blue diamond) and M = 100m (yellow triangle). 
As a guide to the eye, points are connected by solid lines of 
same color. For increasing particle size the KL approaches 
zero, and hence the distribution of the fluctuating force be- 
comes Gaussian. 
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Figure 12: (color online) The deviation of the linear regres- 
sion coefficient from unity, 1 — A is displayed as a function 
of inverse mass 1/M for fluids with densities n* = 0.8 and 
n* = 0.4 in panels (a) and (b), respectively and for Brownian 
particle diameters D = d (black plus), D = l.hd (blue cross), 
D = 2d (red open square), D = 3d (yellow open triangle), 
D = Ad (green open circle), D = 5d (black open diamond), 
D = 7d (blue filled square), and D = lOd (red filled trian- 
gle). For the sake of better readability the data points are 
connected by lines with the color of the corresponding sym- 
bol. In all cases the approach of A to unity for large values 
of mass M is approximately linear. It is faster for smaller 
particles and for the smaller fluid density. 



variance of a distribution p has the minimal diver- 
gence from p compared to all other Gaussian distri- 
butions. The minimal Kullback-Leibler divergence 
of a Gaussian to the fluctuating force distribution 
is displayed in Fig. 11. In particular, it corrobo- 
rates that the fluctuating force indeed approaches 
a Gaussian distribution with increasing diameter of 
the Brownian particle. 

Finally we investigated the relation between the 
total and the fluctuating forces F(t) and F + (t). 
We first discuss the limit of large masses of the 
Brownian particle. Since the memory kernel de- 
creases with increasing mass M whereas the fluctu- 
ating force correlation approaches a finite function 
in this limit, the organized contribution to the to- 
tal force becomes negligible compared to the fluc- 
tuating force. Therefore one expects that the total 
and the fluctuating force agree with each other for 



Brownian particles with sufficiently large mass com- 
pared to the fluid particle mass. In the literature 
[30] one finds as a rough estimate of the deviations 

F+(t) = (1 + 0((m/M)V 2 ))F(t) (33) 

To investigate the relation between F + (t) and F(t) 
in more detail we determined the distributions 
of the fluctuating force conditioned on the total 
force within intervals extending over one force unit. 
Apart from Brownian particles that are small and 
light at the same time it turned out that the distri- 
butions of the conditional fluctuating forces are in- 
dependent of the conditioning value of F(t) except 
for an overall shift which is linear in the condition. 
Therefore the relation between the fluctuating and 
the total force can be described by the model 
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F+ =AF + £ 



(34) 



where Fp denotes the fluctuating force F + con- 
ditioned on the total force restricted to the inter- 
val [F — 0.5e/d, F + Q.he/d). Here, the residue £ is 
independent of the condition F and Gaussian dis- 
tributed. Values of the parameter A were estimated 
by means of linear regression based on eq. (34). De- 
viations of A from unity are displayed in Fig. 12 as 
functions of inverse mass. Apparently A approaches 
unity linearly in 1/M as the mass goes to infinity. 
So actually the convergence of the fluctuating force 
towards the total force is found to be faster than 
indicated by eq. (33). 

8. Conclusions 

We presented an extensive numerical study of 
the motion of single Brownian particles with var- 
ious masses and diameters interacting with fluid 
particles. The data obtained from the MD sim- 
ulations were analyzed in terms of a generalized 
Langevin equation of Mori type. The memory ker- 
nel was estimated on the basis of the total force- 
momentum correlation function and the total force 
auto-correlation functions both of which can di- 
rectly be obtained from the MD simulation data. 
The consistency of the obtained memory kernel and 
the fluctuating forces was tested in terms of the fluc- 
tuation dissipation theorem and an identity for the 
fluctuating force momentum correlation function. 
We discussed the momentum auto-correlation func- 
tion for different Brownian particle masses and di- 
ameters and found pronounced oscillations for large 
and light particles. The maximal value of the mem- 
ory function is inversely proportional to mass. Its 
decay on the macroscopic time scale tm = t<J M/m 
is faster for smaller and/or heavier particles. This 
is in qualitative agreement with Ref. [37] saying 
that the Markovian limit of Brownian dynamics is 
approached if the ratio of fluid and mass densities 
approaches zero. 

The distribution of the fluctuating forces is lep- 
tokurtic for small particles and approaches a Gaus- 
sian form for larger particles. The approach to a 
Gaussian is faster for light particles. For Brownian 
particles being ten times larger than the fluid parti- 
cles the fluctuating force was found to be perfectly 
Gaussian for all considered masses up to 100 times 
the fluid particle mass. This observation is based 
on direct comparison of histograms, on estimates of 
the kurtosis and on the Kullback-Leibler divergence 
of the closest Gaussian. 



For massive Brownian particles the difference be- 
tween the fluctuating and the total force was shown 
to shrink proportionally to the mass ratio m/M. 
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A. Proof of eq. (29) and the fluctuation dis- 
sipation theorem 

We first show that eq. (28) in combination with 
eq. (17) implies eq. (29). Solving the stationary 
generalized Langevin equation (28) for the fluctu- 
ating force yields 

F + (t)=F(t)+ I ds k(t - s)P(s) , (35) 

J — OQ 

where the total force F(t) gives rise to the acceler- 
ation of the Brownian particle, i.e. 

P(t) = F(t) . (36) 

By multiplying both sides of eq. (35) with P(0) and 
performing an equilibrium average we obtain 

(F+(t)P) = (P(t)P(O)) + [ ds k(t-s)(P(s)P(0)) 

J — DC 

= C(t)+[ dsk(t-s)C(s)+ [ dsk(t-s)C(s) . 

JO J-oc 

(37) 

The first two terms of the right hand side of the last 
line cancel each other as a consequence of eq. (17) 
such that only the third term remains. Changing in 
this remaining term the variable of integration s — > 
— s and observing the symmetry of the correlation 
function (P(t)P(O)) = C(t) = C(-t) one obtains 
the expression on the right hand side of eq. (29). 
This proves eq. (29). 



In order to prove that the fluctuation dissipa- 
tion theorem follows from the generalized Langevin 
equation (28) in combination with the form of the 
fluctuating force-momentum correlation cq. (29) 
and the equation of motion of the momentum auto- 
correlation function eq. (17) we determine the auto- 
correlation function of the fluctuating force using 
eq. (35). We then obtain with eq. (35) 

(F+(t)F+(0)) = (F+(t)F(Q)) 



dsk(-s)(F+(t)P(s)) 



-c(t) 



ds k(t-s)C(s) 



r o ,oo (38) 
+ / ds ds' k(-s)k(t-s+s')C(s') 



-oo JO 



fc(t)C(0)+ / dsk(t + s)C(s) 
Jo 



oo c oo 



ds 



ds' k(s)k(t 



s')C(s') 



where we used eq. (29) at the second equal sign 
and eq. (19) at the third equal sign. Using eq. (17) 
one may show with little algebra that the last two 
terms compensate each other such that the fluctu- 
ation dissipation theorem is recovered. 

Vice versa, if one takes for granted the stationary 
generalized Langevin equation (28) together with 
the fluctuation dissipation theorem (8) and the fluc- 
tuating force- momentum correlation function (29), 
one finds for the momentum auto-correlation func- 
tion the following equation of motion 



C(t) 



— oo 

i 



dsk(t-s)C(s) + (F+(t)P(0)) 
ds k(t - s)C(s) , 



(39) 

where we split the integral on the right hand side 
of the first line into a part extending from — oo to 
and another one from to t. The first contri- 
bution is exactly canceled by the fluctuating force- 
momentum correlation given by eq. (29). Hence we 
recover eq. (17) as the equation of motion for the 
momentum auto-correlation function. 
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